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Surfaces eroded by ion-sputtering are sometimes observed to develop morphologies which are 
either ripple (periodic), or rough (non-periodic). We introduce a discrete stochastic model that 
allows us to interpret these experimental observations within a unified framework. We find that 
a periodic ripple morphology characterizes the initial stages of the evolution, whereas the surface 
displays self-afnne scaling in the later time regime. Further, we argue that the stochastic continuum 
equation describing the surface height is a noisy version of the Kuramoto-Sivashinsky equation. 
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A remarkable feature of erosion processes via ion sput- 
tering in amorphous materials is the formation of a pat- 
tern consisting of a ripple structure, aligned in directions 
either parallel to or perpendicular to that of the bom- 
barding beam of ions |,|. Indeed, we might expect that 
erosion tends to erase every possible feature of the sur- 
face morphology, and that the presence of noise in the 
system would further act against the formation of such a 
periodic pattern. Only recently have there been experi- 
mental H and theoretical Q attempts to understand the 
formation of a ripple structure in the more general con- 
text of non-equilibrium interface growth phenomena j| . 
Many such interfaces are "rough" and exhibit self-affinc 
scaling at long distances and long times M; experimen- 
tally, one finds that surfaces eroded by ion bombardment 
also exhibit self-affine scaling behavior || . An outstand- 
ing question is then how to reconcile these observations 
with the formation of the periodic ripple structure. 

In this Letter, we introduce a discrete stochastic model 
that incorporates the main physical mechanisms believed 
to influence the dynamics of the eroded surface morphol- 
ogy. We show that this model is characterized by an 
initial stage in which the surface morphology displays a 
ripple structure, and that subsequent stages are charac- 
terized by a crossover to a rough surface in the univer- 
sality class of the Kardar-Parisi-Zhang (KPZ) equation 
||. We argue that the stochastic equation which pro- 
vides the continuum description of our model is a noisy 
version of the Kuramoto-Sivashinsky (KS) equation 
Thus we interpret the formation of a periodic pat- 
tern H and the development of rough interfaces [£| in 
the ion-sputtered systems as the early and late regimes 
respectively of the same dynamical process. 

We first consider the main mechanisms that deter- 
mine the surface morphology undergoing ion bombard- 
ment. 

(i) Erosion — In ion sputtering, an initially flat sub- 
strate is bombarded with a well-collimated beam of heavy 
ions carrying a certain kinetic energy, and forming a pre- 
cise angle with the normal to the uneroded surface. The 
phenomena leading to erosion take place within some fi- 



nite distance from the surface. Namely, the ions pene- 
trate inside the solid and induce along their path cas- 
cades of collisions among the atoms of the substrate. 
Atoms located at the surface may be affected by these 
collisions and acquire enough energy to leave the sur- 
face Consequently, atoms located at the bottom of 
troughs gain more energy on average and they are pref- 
erentially eroded as compared to those on the peaks of 
crests |jl^] . This "instability" can be thought of as a neg- 
ative surface tension, since the surface tends to maximize 
its area. 

(ii) Surface Diffusion — In physical systems there is 
a stabilizing mechanism that balances the negative sur- 
face tension, surface diffusion, which is always present at 
a non-zero temperature. Particles on the surface tend 
to diffuse looking for highly coordinated positions, a 
relevant phenomenon for sputtering |5]|^] as well as for 
molecular-beam epitaxy (MBE) ( [l3| . 

When the angle between the local normal to the sur- 
face and the incident beam approaches the grazing value, 
there is an increase in the reflection of the ions by the 
surface and the rate of erosion diminishes. This surface 
effect is beyond the approximations made in Sigmund's 
theory of ion sputtering O] , and is reflected in the angle 
dependence of the sputtering yield, Y(tp), defined to be 
the number of eroded particles divided by the total num- 
ber of bombarding ions. Here tp is the local angle of the 
ion trajectories to the surface normal at each point p4[ . 
Typically Y(ip) is symmetric around <p = 0, presents 
a maximum between 60° and 80°, and decreases to zero 
as ip -> 90°. 

To define the model, we introduce two dynamical rules, 
one to account for erosion and one to account for surface 
diffusion. The rule for erosion incorporates the unstable 
behavior described above as well as the phenomenolog- 
ical dependence of the sputtering yield as a function of 
ip. The model for the case of 1 + 1 dimensions |1| is 
defined on a square lattice of lateral size L, with periodic 
boundary conditions in the horizontal direction. The ini- 
tial interface is a horizontal line separating occupied sites 
(below) from empty sites (above). We choose randomly a 
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site % at the interface where % = 1, . . . , L. The chosen site 
is subject to erosion with probability p, or to diffusion 
with probability 1 — p, where the rules are as follows: 

(i) Erosion (probability p) — Jl(f | We compute tp = 
tan _1 [(/i,; + i — hi^i)/2], where hi is the height of the inter- 
face at site i, and apply the erosion rule with probability 
Y(ip), as given in Fig. la. To erode, we count the number 
of occupied neighbors inside a square box of size 3x3 
lattice spacings centered in the chosen site i (box rule). 
We empty the site with an erosion probability p e pro- 
portional to the number of occupied cells in the box (see 
Fig. 16). Thus the box rule favors the erosion of troughs 
as compared to the peaks of crests, and therefore is the 
source of the instability in the ion-sputtered system. 

(ii) Surface Diffusion (probability 1 — p) — A diffu- 
sive move of the particle i to a nearest neighbor col- 
umn is attempted with hopping probability Wi—>f = 
[1 + exp (AHi^f /kgT)] -1 , where ATL^f is the energy 
difference between the final and initial states of the move. 
Following 0, we choose H = (J/2) Y,(i,j) {hi ~ h ) 2 . 

A continuum equation that describes the dynamics of 
the interface height has the form 

d t h(x,t) = vV 2 h- nV 4 h + n{x,t) + f Y [h(x,t)}. (1) 

Here h(x,t) is the height of the interface at position x 
and time t, v is a negative surface tension coefficient, k 
is a positive coefficient that accounts for the surface dif- 
fusion, and rj(x,t) is a Gaussian noise term with short 
range correlations and strength 2D, that accounts for 
the fluctuations in the flux of incoming particles. The 
functional fy[h] takes into account the contribution of 
nonlinear terms, which appear in the equation of motion 
due to the effect of Y ((f), itself a nonlinear function of the 
local slope V/i = tan if. The nonlinearities become more 
relevant in the equation of motion at the late regime of 
the evolution when large slopes develop, see below. 

First we study the model in the case Y(ip) = 1, 
which approximately holds at the early stages of evo- 
lution. We focus on the time dependence of the total 
interface width W{t) = {L' 1 £f =1 ~ W)Y) 1,2 i 

where hit) = L^ 1 ^^ =i^i(*)> an( i the brackets denote 
an average over realizations of the noise. The erosion 
rule alone — corresponding to p — 1, Y(ip) = 1 — leads 
to W(t) ~ t, which can be attributed to v in ([j]) being 
a negative number pgj |. Moreover, by considering only 
the erosion rule, but defined with probability 1 — p e , we 
find that the interface has the scaling properties of the 

that is obtained 
h] = 0. We un- 



Edwards- Wilkinson (EW) equation [19 
from (0) with v > 0, k = 0, and fy 
derstand this result since one favors the erosion of peaks 
as compared to valleys, leading to the smoothing mecha- 
nism characteristic of a positive surface tension. On the 
other hand, the surface diffusion mechanism alone (for 
p = 0, Y(ip) = 1) is well described by (^) with v = 0, 
k > 0, and fy[h] = 0, i.e., the linear MBE equation 
jl7],[l3]]. When mechanisms (i) and (ii) above are con- 
sidered simultaneously for J/ksT = 5, p = 0.5, and 



Y((p) = 1, we obtain the different stages of the time evo- 
lution displayed in Fig. 2a |^lj . There exists a first region 
H for which W{t) ~ tP\ with /3i = 0.38 ± 0.02, the 
growth exponent for the linear MBE universality class, 
after which W(t) scales with $2 > 0.5 due to the insta- 
bility caused by v < 0. In this linear stability 
analysis of ([j]) shows that there is a maximally unstable 
mode in the system, k m = (\v\/2k) 1 ^ 2 , and therefore the 
surface is almost periodic (inset in Fig. 2a). 

Next we consider the model with Y((p) shown in Fig. 
la. The results are not expected to depend strongly on 
the specific form of Y(ip), so long as it preserves the exis- 
tence of a maximum, and Y(0) / 0, Y(90°) = Fig. 

shows the structure factor S(k) = (h(k,t)h(—k,t)) at 
the onset of the instability. Here h(k, t) is the Fourier 
transform of hi(t) — h(t). As we see, the early stages 
of the dynamics are still dominated by the periodic rip- 
ple structure defined by the competition between surface 
tension and surface diffusion, described by the linear part 
of®. 

For later times, the large slopes built up by the in- 
stability induce nonlinear effects, and the interface re- 
sults in a rough morphology (inset of Fig. 2b). In Fig. 
2b, we present the time evolution of W(t) for the com- 
plete model. We again observe a first regime [^2j with 
0i = 0.38 ±0.03, followed by unstable erosion (/3 2 > 0.5). 
For later stages, we find = 0.23 ± 0.03, consistent 
with EW, after which a crossover to (3^ = 0.28 ± 0.03 
is found. Finally, the width saturates due to the finite 
size of the system. Note that the value of the growth 
exponent for the KPZ equation is /3kpz = 1/3 0. At 
saturation, S(k) displays the small momenta behavior 
S{k, t) ~ fc -2 , consistent with the scaling of both the EW 
and KPZ universality classes (see Fig. 3). To determine 
if a KPZ nonlinearity is present in Eq. (0), we compute 
the mean velocity v(m) of the interface in the saturated 
regime as a function of an average tilt m = (V/i) im- 
posed by using helical boundary conditions. If we as- 
sume that the relevant nonlinearity in (Q) is of the KPZ 
type, then fy\h] — (A/2)(V/i) 2 . Taking spatial and noise 
averages in ([j]), v = vq + (A/2)m 2 , where vq is the ve- 
locity of the untilted interface |24| . The parabolic shape 
of v(m) obtained in our simulations (see Fig. 4) leads to 
the conclusion that the long time and long distance be- 
havior of the model falls into the KPZ universality class. 
Moreover, the continuum equation describing the model 
ion-sputtered surfaces is the noisy KS equation 



dth= vV 2 h~n\J A h 



^(Vh) 2 + Tj(x,t). (2) 



To compare the dynamics of (g) with those obtained 
for the discrete model, we integrate numerically Eq. (||) 
in 1 + 1 dimensions. Fig. || shows the behavior of the 
function W{t). We observe the same crossovers as in the 
model |2^] : f3\ — 0.38 ± 0.03 corresponding to the lin- 
ear MBE case, followed by unstable growth (ii > 0.5. 
Then, a transition to EW behavior f3 3 = 0.25 ± 0.03 
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is observed, after which the nonlinearities dominate and 
a KPZ growth with /3 4 = 0.30 ± 0.03 is obtained. Fi- 
nally, the interface width saturates due to the finite size 
of the system. Consistent with these numerical findings, 
the late scaling of Eq. (|J) has been shown through a 
renormalization-group calculation to be that of the 
KPZ equation in 1+1 and 2+1 dimensions. As we see, 
both in the model and the noisy KS, there is a long 
crossover time from EW to KPZ behavior, responsible 
for the difference between /?4 and /3k pz, and for the nar- 
row window in which (3^ is observed — we find that the 
width of this window increases systematically with L. A 
similar phenomenon is well known to occur in the deter- 
ministic KS equation in 1+1 dimensions, see Sneppen et 
al. in [§. 

Finally, we compare the results of the model with ob- 
servations of recent experiments pal . The experimental 
development of a ripple structure [f| is well understood 
in terms of the unstable linear theory of ion-sputtering 
describing the early stages of the time evolution of the 
model presented here. Moreover, the model predicts that 
in the late regime the large slopes generated by the un- 
stable growth trigger the action of nonlinearities which 
stabilize the surface. The nonlinearity we find is of the 
KPZ type, consistent with the experimental observation 
of KPZ scaling reported by Eklund et al. ||. To confirm 
the above picture, it would be of interest to study exper- 
imentally if both regimes do effectively take place in the 
time evolution of the same physical system. 
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FIG. 1. (a) Sputtering yield Y(ip) as a function of the an- 
gle cp. (b) Box rule for erosion. We define p e as the number 
of occupied neighboring sites (grey squares) inside the 3x3 
box centered at site i (black square), normalized by 7. The 
examples shown correspond to (i) p e = 1 and (ii) p e = 3/7. 

FIG. 2. (a) Time evolution of the surface width for the 
cases Y(ifi) = 1 and L — 50. The solid line is the consecutive 
slope of the width, showing the value of the growth exponent 
/3 in each regime. The inset shows the ripple structure of the 
interface at t — 1000. The saturation observed in W(t) is 
due to the discreteness of the lattice: the erosion rule breaks 
down when the local slopes of the interface are bigger than 
3 [^. This effect can be avoided by using a bigger box. 
(6) Interface width as a function of time for the full model 
showing the regimes of the evolution for L = 2048. As in 
(a), the solid line is the consecutive slope. The inset shows 
a portion of the rough morphology at the late regime, where 
the self-afnne scaling behavior holds. The arrows indicate the 
times at which the structure factor is displayed in Fig. 3. 



FIG. 3. Structure factor computed using the full model for 
a system with L = 2048. For t = 300, averaged over 2600 
noise realizations (o), and for t = 1.7 x 10 6 , averaged over 39 
realizations (•) (see arrows in Fig. 2b). The solid line is a fit 
to the exact solution of the discretized linear part of Eq. (1). 
The dashed straight line has slope —2. 

FIG. 4. Plot of the mean velocity v(m) as a function of 
the average tilt m of the interface, calculated in the saturated 
regime for L — 128, and averaged over 810 noise realizations. 
The solid line is a fit to a parabola. 

FIG. 5. Log-log plot of the interface width as a function of 
time obtained in the numerical integration of Eq. (2). The 
mesh is Ax = 1.0 and the time step At = 0.005. Parameter 
values are v — — 1, k = 1, A = 4, D = 0.5, and L — 1024. 
The solid line gives the consecutive slopes. 



4 




1 0" 2 1 0" 1 10 



k 



-0.21425 r 
-0.21460 - 
-0.21495 - 
-0.21530 - 
-0.21565 - 
-0.21600 - 



